Le leader mondial ‘Talend’ de l’intégration et de l’intégrité des données dans le cloud, a annoncé le 30 juillet 2019 la disponibilité de Stitch Data Loader sur la marketplace d’Amazon Web Services (AWS). Déjà pertinente par sa portée et rendue d’autant plus accessible depuis cette annonce, la plateforme Stitch est devenue un incontournable dans la procuration de données ouvertes. Stitch Data Loader facilite en effet l’ingestion de données pour l’analyse et s’adapte automatiquement à l’augmentation du volume de données.
Cette plateforme est une composante de Talend Data Fabric, qui est une suite d’applications cloud conçues pour aider les clients à collecter, gérer, transformer et partager données et informations.
Pour ce projet, ont été récupérées différentes données issues de la plateforme Stitch concernant différentes requêtes effectuées, appelées « jobs » dont nous chercherons à prédire la durée en secondes.
Le but de ce projet est de créer des modèles additifs généralisés qui prédisent au mieux le temps de durée de chaque job sur la plateforme Stitch. Nous construirons pour cela différents GAM sur la base d’entraînement, modèles que nous chercherons à optimiser pour les rendre les plus performants possible, puis nous les appliquerons sur la base ‘test’ où nous espérons les erreurs les plus faibles.
Dans un premier temps nous appelons les librairies qu’il nous faut pour réaliser le projet :
library(tidyverse)
library(readr)
library(gganimate)
library(viridis)
library(gridExtra)
library(summarytools)
library(kableExtra)
library(knitr)
library(flextable)
library(tsoutliers)
library(EnvStats)
library(stringr)
library(seastests)
#install.packages("rAmCharts")
library(rAmCharts)
library(mgcv)
library(car)
Nous avons téléchargé les données depuis la plateforme Kaggle :
train <- read_csv('./Données/data_train.csv')
test <- read_csv('./Données/data_test.csv')
Comme évoqué précédemment, les données brutes contiennent 5 variables et 9100 observations.
Nous disposons initialement de 6 variables :
Nous regardons les 5 premières lignes de notre jeu de données pour nous rendre compte des informations qu’il contient. Nous pouvons aussi regarder le format des variables grâce à la commande str(). Puis, le summary de la base nous donnera un aperçu des distributions des variables avec les valeurs minimales et maximales de chacune d’elle ainsi que la moyenne, les quartiles etc. Le but de cette analyse exploratoire étant de cerner au mieux les données pour éventuellement créer de nouvelles variables ou transformer les existantes.
| CONNECTION_ID | START_TIME | TAP_VERSION | TAP_EXIT_STATUS | JOB_DURATION |
|---|---|---|---|---|
| 1 | 2019-10-11 13:41:00 | 1.15.1 | 0 | 9 |
| 1 | 2019-10-11 14:41:00 | 1.15.1 | 0 | 9 |
| 1 | 2019-10-11 15:41:00 | 1.15.1 | 0 | 12 |
| 1 | 2019-10-11 16:41:00 | 1.15.1 | 0 | 10 |
| 1 | 2019-10-11 17:41:00 | 1.15.1 | 0 | 11 |
| 1 | 2019-10-11 18:41:00 | 1.15.1 | 0 | 10 |
| 1 | 2019-10-11 19:41:00 | 1.15.1 | 0 | 10 |
| 1 | 2019-10-11 20:41:00 | 1.15.1 | 0 | 11 |
| 1 | 2019-10-11 21:41:00 | 1.15.1 | 0 | 10 |
| 1 | 2019-10-11 22:41:00 | 1.15.1 | 0 | 11 |
On constate d’emblée que notre variable à expliquer “JOB_DURATION” correspond à la cinquième colonne, et que la variable “START_TIME” regroupe à la fois les dates et les heures du début des jobs. À partir des 10 premières observations, on voit que les variables “CONNECTION_ID”, “TAP_VERSION” et “TAP_EXIT_STATUS” restent inchangées, ce qui n’est pas le cas de “JOB_DURATION” qui évolue en fonction du temps récupéré par la variable “START_TIME”. De plus, les durées des 10 premiers jobs sont très courtes, allant de 9 à 12 secondes.
str(train)
## Classes 'spec_tbl_df', 'tbl_df', 'tbl' and 'data.frame': 9100 obs. of 5 variables:
## $ CONNECTION_ID : num 1 1 1 1 1 1 1 1 1 1 ...
## $ START_TIME : POSIXct, format: "2019-10-11 13:41:00" "2019-10-11 14:41:00" ...
## $ TAP_VERSION : chr "1.15.1" "1.15.1" "1.15.1" "1.15.1" ...
## $ TAP_EXIT_STATUS: num 0 0 0 0 0 0 0 0 0 0 ...
## $ JOB_DURATION : num 9 9 12 10 11 10 10 11 10 11 ...
## - attr(*, "spec")=
## .. cols(
## .. CONNECTION_ID = col_double(),
## .. START_TIME = col_datetime(format = ""),
## .. TAP_VERSION = col_character(),
## .. TAP_EXIT_STATUS = col_double(),
## .. JOB_DURATION = col_double()
## .. )
Sur nos 5 variables seules 3 sont numériques ; “TAP_VERSION” est en format ‘character’ et “START_TIME” au format ‘POSIXct’ ce que nous devrons donc changer pour un meilleur traitement des données.
summary(train)
## CONNECTION_ID START_TIME TAP_VERSION
## Min. :1.000 Min. :2019-10-11 13:41:00 Length:9100
## 1st Qu.:1.000 1st Qu.:2020-01-11 05:31:00 Class :character
## Median :1.000 Median :2020-02-27 21:50:30 Mode :character
## Mean :1.626 Mean :2020-03-19 21:41:01
## 3rd Qu.:2.000 3rd Qu.:2020-07-13 01:46:30
## Max. :3.000 Max. :2020-08-29 11:09:00
##
## TAP_EXIT_STATUS JOB_DURATION
## Min. :0.00000 Min. : 2.0
## 1st Qu.:0.00000 1st Qu.: 10.0
## Median :0.00000 Median : 14.0
## Mean :0.00099 Mean : 400.5
## 3rd Qu.:0.00000 3rd Qu.: 181.0
## Max. :1.00000 Max. :4349.0
## NA's :11
On voit que sur les 3 modalités de connexion, la 1 regroupe au moins la moitié des observations puisque la médiane est de 1. On constate aussi que les requêtes de la base train s’étendent du 11 octobre 2019 au 29 août 2020. Concernant la variable “TAP_EXIT_STATUS” c’est la seule qui contient des valeurs manquantes, il faudra donc retirer ces 11 observations pour les modélisations. Enfin, on voit que les durées des requêtes vont de 2 secondes à 4349 soit 1 heure 12 minutes et 29 secondes, ce qui est considérable. On peut donc supposer que le job n’a pas abouti (TAP_EXIT_STATUS=1) ou bien qu’il existe une grande disparité des temps au sein des périodes. On voit ainsi que la moyenne est de 400 secondes tandis que la médiane est de seulement 14.
Avant toute transformation des variables nous allons visualiser la variable à expliquer en fonction du temps.
# On plot la variable à expliquer fonction du temps
plot(train$JOB_DURATION ~ train$START_TIME, type='l', xlab="Temps", ylab="Durées des jobs en secondes", main="Évolution de la durée des requêtes dans le temps")
Au vu du graphique ci-dessus on peut constater qu’il y a 3 périodes correspondant aux 3 chroniques :
La connexion ID n°3 se distingue donc largement des chroniques 1 et 2 par sa grande volatilité de temps des jobs, malgré le fait qu’elle ait été capturée sur une sous-période de la chronique 1. Compte tenu du fait que les valeurs des jobs diffèrent beaucoup d’une chronique à l’autre, nous allons procéder à la division de notre base en 3 sous bases pour chaque chronique, dans le but d’observer plus précisément les variations des temps des jobs.
# On créé une sous-base pour la chronique n°1
df_ID1 <- train %>% filter(CONNECTION_ID==1)
summary(df_ID1)
## CONNECTION_ID START_TIME TAP_VERSION TAP_EXIT_STATUS
## Min. :1 Min. :2019-10-11 13:41:00 Length:5000 Min. :0
## 1st Qu.:1 1st Qu.:2019-12-02 14:26:00 Class :character 1st Qu.:0
## Median :1 Median :2020-01-23 16:11:00 Mode :character Median :0
## Mean :1 Mean :2020-01-23 22:53:32 Mean :0
## 3rd Qu.:1 3rd Qu.:2020-03-16 09:56:00 3rd Qu.:0
## Max. :1 Max. :2020-05-07 12:41:00 Max. :0
## JOB_DURATION
## Min. : 7.00
## 1st Qu.: 10.00
## Median : 10.00
## Mean : 12.25
## 3rd Qu.: 12.00
## Max. :128.00
La sous base de la chronique 1 ainsi créée contient 55% de la base train pour un total de 5000 observations. Nous constatons que pour la variable “TAP_EXIT_STATUS” il n’y a plus de valeur manquante, celles-ci devaient donc concerner les chroniques 2 ou 3 (nous verrons cela juste après). En revanche si l’on regarde la durée des jobs, on voit que le temps maximal est de 128 secondes, il s’agit sûrement d’une valeur atypique qui peut gêner nos modélisations et que nous allons donc traiter dès maintenant.
amBoxplot(df_ID1$JOB_DURATION, xlab="CONNECTION ID n°1", ylab="Durée des jobs", main="Boxplot 'JOB_DURATION' pour la chronique n°1")
À partir du graphique ci-dessus, nous observons d’emblée que de nombreuses observations sont supérieures au troisième quantile et peuvent donc être considérées comme atypiques. Cependant, étant donné le nombre important de points concernés il peut s’agir d’un comportement explicable et qui peut apporter une information importante pour les futures modélisations. Nous considérons donc comme atypiques les valeurs qui dépassent 1 minute car il n’y en a que 8 qui sont assez espacées, et procédons au test de Rosner pour confirmer leur atypicité.
rosnerTest(df_ID1$JOB_DURATION, k = 8, alpha = 0.05) # 8 outliers : 128,126,121,114,108,76,75,65
# On regarde à quelles dates les points correspondent pour voir si pattern
a <- df_ID1 %>% slice(1423,846,967,2479,487,1711,366,2489)
kbl(a, align = "c") %>% kable_classic(full_width = F, html_font = "Cambria")
| CONNECTION_ID | START_TIME | TAP_VERSION | TAP_EXIT_STATUS | JOB_DURATION |
|---|---|---|---|---|
| 1 | 2019-12-09 18:41:00 | 1.15.1 | 0 | 128 |
| 1 | 2019-11-15 17:41:00 | 1.15.1 | 0 | 126 |
| 1 | 2019-11-20 18:41:00 | 1.15.1 | 0 | 121 |
| 1 | 2020-01-22 18:41:00 | 1.15.1 | 0 | 114 |
| 1 | 2019-10-31 18:41:00 | 1.15.1 | 0 | 108 |
| 1 | 2019-12-21 18:41:00 | 1.15.1 | 0 | 76 |
| 1 | 2019-10-26 18:41:00 | 1.15.1 | 0 | 75 |
| 1 | 2020-01-23 04:41:00 | 1.15.1 | 0 | 65 |
D’après le test de Rosner on voit que les 8 observations étaient effectivement atypiques, nous trouvions intéressant de regarder de quelles observations il s’agissait pour noter d’éventuels comportements communs par exemple par heure ou par jour de la semaine. Nous observons alors d’après la table que la plupart des observations atypiques se situent en fin d’après-midi à 18h41 (c’est le cas de 6 observations sur 8). On ne note cependant pas de régularité dans les jours de la semaine puisque l’on retrouve 2 fois mercredi, jeudi et samedi, et 1 fois lundi et vendredi.
# on crée une nouvelle base sans ces points
df_ID1=df_ID1[-c(1423,846,967,2479,487,1711,366,2489),]
dim(df_ID1)
## [1] 4992 5
amBoxplot(df_ID1$JOB_DURATION, xlab="CONNECTION ID n°1", ylab="Durée des jobs", main="Boxplot 'JOB_DURATION' débarrassé des points atypiques")
On retire donc les 8 observations de la base df_ID1, et en effectuant à nouveau un diagramme à moustaches on observe que les observations sont assez rapprochées les unes des autres et peuvent ainsi nous apporter de l’information (en cliquant dessus on voit que les valeurs les plus élevées hors de la boîte sont à 58 secondes). Nous allons à présent regarder les temps de durée des requêtes Stitch pour la chronique 1 en visualisant pour toute la base, par journée, par semaine et par mois; l’objectif étant de noter une possible saisonnalité pour la variable à expliquer “JOB_DURATION”.
# Graphiques CONNECTION_ID 1
par(mfrow=c(2,2))
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1, main="Toute la période", xlab="Temps", ylab="Durées des jobs") # le tout
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[84:107,], main="Un mardi", xlab="Temps", ylab="Durées des jobs") # mardi
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[564:731,], main="Une semaine", xlab="Temps", ylab="Durées des jobs") # une semaine
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[491:1208,], main="Un mois (novembre)", xlab="Temps", ylab="Durées des jobs") # un mois
On peut constater à partir des 4 graphiques des tendances :
par(mfrow=c(2,2))
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[60:227,], main="Une semaine", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[228:395,], main="Une semaine", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[396:563,], main="Une semaine", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[564:731,], main="Une semaine", xlab="Temps", ylab="Durées des jobs")
par(mfrow=c(2,2))
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[491:1208,], main="Un mois (novembre)", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[1209:1926,], main="Un mois (décembre)", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[1927:2644,], main="Un mois (janvier)", xlab="Temps", ylab="Durées des jobs")
plot(JOB_DURATION ~ START_TIME, type='l', data=df_ID1[2645:3362,], main="Un mois (février)", xlab="Temps", ylab="Durées des jobs")
# On créé une sous-base pour la chronique n°2
df_ID2 <- train %>% filter(CONNECTION_ID==2)
summary(df_ID2)
## CONNECTION_ID START_TIME TAP_VERSION
## Min. :2 Min. :2020-07-08 09:39:00 Length:2500
## 1st Qu.:2 1st Qu.:2020-07-21 10:01:30 Class :character
## Median :2 Median :2020-08-03 10:24:00 Mode :character
## Mean :2 Mean :2020-08-03 10:24:00
## 3rd Qu.:2 3rd Qu.:2020-08-16 10:46:30
## Max. :2 Max. :2020-08-29 11:09:00
##
## TAP_EXIT_STATUS JOB_DURATION
## Min. :0.000000 Min. : 67.0
## 1st Qu.:0.000000 1st Qu.:165.8
## Median :0.000000 Median :172.0
## Mean :0.001601 Mean :178.4
## 3rd Qu.:0.000000 3rd Qu.:182.0
## Max. :1.000000 Max. :688.0
## NA's :1
Nous réitérons l’analyse pour la chronique 2 en créant une autre sous base pour cet ID qui contient cette fois 2500 observations. On voit cette fois que la durée moyenne des jobs est beaucoup plus importante que pour la chronique 1 variant autour de 178 secondes. De plus nous observons qu’il y a une valeur manquante pour la variable “TAP_EXIT_STATUS” que nous retirons avant la détection des outliers.
#on supprime la valeur manquante pour "TAP_EXIT_STATUS" (6 aout 2020)
df_ID2 <- na.omit(df_ID2)
#on regarde le boxplot
amBoxplot(df_ID2$JOB_DURATION, xlab="CONNECTION ID n°2", ylab="Durée des jobs", main="Boxplot 'JOB_DURATION' pour la chronique n°2")
| CONNECTION_ID | START_TIME | TAP_VERSION | TAP_EXIT_STATUS | JOB_DURATION |
|---|---|---|---|---|
| 2 | 2020-07-30 12:09:00 | 1.4.33 | 0 | 688 |
| 2 | 2020-07-28 17:39:00 | 1.4.33 | 0 | 604 |
Comme visible sur le boxplot il y a 2 valeurs extrêmes correspondant aux observations 1062 et 977 (détectées avec le test de Rosner) ; les 2 durées atypiques d’environ 10 minutes d’attente ont eu lieu à deux journées d’intervalle fin juillet 2019. Les temps de la variable “JOB_DURATION” sont maintenant compris entre 151 et 396 secondes. En moyenne le temps d’attente pour la chronique 2 est 15 fois plus long que pour la chronique 1.
Par rapport aux graphiques de la chronique une, on retrouve le pic du temps des jobs de début de journée assez grand et un pic plus faible en fin d’après-midi. En revanche on constate que le pic de début de journée est légèrement plus tard qu’en période “scolaire” puisque nous le rappelons la chronique 2 concerne les mois de juillet et août 2020.
# On créé une sous-base pour la chronique n°3
df_ID3 <- train %>% filter(CONNECTION_ID==3)
summary(df_ID3)
## CONNECTION_ID START_TIME TAP_VERSION
## Min. :3 Min. :2020-01-08 03:02:00 Length:1600
## 1st Qu.:3 1st Qu.:2020-01-24 17:45:00 Class :character
## Median :3 Median :2020-02-10 09:30:00 Mode :character
## Mean :3 Mean :2020-02-10 10:02:12
## 3rd Qu.:3 3rd Qu.:2020-02-27 01:15:00
## Max. :3 Max. :2020-03-14 19:00:00
##
## TAP_EXIT_STATUS JOB_DURATION
## Min. :0.000000 Min. : 2
## 1st Qu.:0.000000 1st Qu.:1664
## Median :0.000000 Median :1824
## Mean :0.003145 Mean :1961
## 3rd Qu.:0.000000 3rd Qu.:2165
## Max. :1.000000 Max. :4349
## NA's :10
Finalement, nous réalisons une troisième fois l’analyse pour la dernière connexion : la base contient alors 1600 observations et 10 valeurs manquantes que nous enlevons de la base. Pour cette période de janvier à mars 2020 on voit que les durées des requêtes sont encore plus élevées que la chronique 2.
| CONNECTION_ID | START_TIME | TAP_VERSION | TAP_EXIT_STATUS | JOB_DURATION |
|---|---|---|---|---|
| 3 | 2020-03-05 21:00:00 | 2.6.5 | 0 | 4349 |
| 3 | 2020-01-08 06:01:00 | 2.6.4 | 0 | 660 |
| 3 | 2020-01-16 05:01:00 | 2.6.4 | 1 | 841 |
| 3 | 2020-01-21 20:00:00 | 2.6.4 | 0 | 562 |
| 3 | 2020-02-09 07:01:00 | 2.6.4 | 1 | 16 |
| 3 | 2020-03-02 01:01:00 | 2.6.5 | 1 | 7 |
| 3 | 2020-03-13 01:01:00 | 2.6.5 | 1 | 8 |
À partir du boxplot on voit que les observations qui semblent atypiques sont celles qui ont un temps d’éxécution inférieur à 1000 et supérieur à 3000 secondes, nous décidons pour notre part de retirer les valeurs non comprises dans l’intervalle [1000;4000] où 7 valeurs sont concernées. Il y a ainsi 3 points en janvier, 1 en février et 3 autres en mars - valeurs que nous retirons de la base df_ID3. En moyenne sur cette sous-période de la chronique une, les jobs durent 1969 secondes, soit plus de 30 minutes.
On peut constater plusieurs choses cette fois-ci :
Cette partie nous a permis de constater la proximité des temps de jobs pour les chroniques 1 et 2 quoiqu’un peu plus longues pour la chronique 2 du fait peut être de l’été. En revanche les temps d’attente des jobs de la chronique 3 sont eux beaucoup plus importants avec des pics à des horaires différents des autres CONNECTION_ID.
Nous allons procéder à la transformation des variables qui nous ont semblées pertinentes dans l’explication du temps des jobs. Ayant découpé l’échantillon ‘train’ en 3 sous-bases selon la chronique, la seule variable qui reste pertinente pour l’explication des jobs est la variable “START_TIME”. A partir de cette dernière nous allons pouvoir extraire plusieurs informations pour expliquer le temps des jobs. En découpant “START_TIME” en nombre d’heures par jour, des “paterns” ont été mis en évidence quelles que soit la chronique. Lorsque nous avons découpé cette variable par semaine pour voir si le jour de la semaine pouvait influencer le temps des jobs, nous nous sommes rendu compte que seule la chronique 3 semblait être affectée par le samedi et dimanche que nous avons appelé précédemment “effet week-end”. Enfin, lorsque l’on a observé Y sur plusieurs mois nous n’avons pas noté de ressemblances ou de tendances d’un mois à l’autre. Dans tous les cas nous faisons le choix d’extraire de la variable “START_TIME” l’heure de la journée qui nous semble être la variable la plus importante, le jour de la semaine et la semaine de l’année, tels que :
# On crée toutes ces variables en séparant chaque info de 'Start_time'
# Semaine de l'année
df_ID1$Week <- format(df_ID1$START_TIME, format = "%U")
df_ID2$Week <- format(df_ID2$START_TIME, format = "%U")
df_ID3$Week <- format(df_ID3$START_TIME, format = "%U")
# Jours de la semaine
#on récupère le nom de la semaine
df_ID1$Day <- format(df_ID1$START_TIME, format = "%A")
df_ID2$Day <- format(df_ID2$START_TIME, format = "%A")
df_ID3$Day <- format(df_ID3$START_TIME, format = "%A")
#on convertit le nom en chiffre
df_ID1$Day <- str_replace_all(df_ID1$Day, c("Lundi"="1", "Mardi"="2", "Mercredi"="3", "Jeudi"="4", "Vendredi"="5", "Samedi"="6", "Dimanche"="7"))
df_ID2$Day <- str_replace_all(df_ID2$Day, c("Lundi"="1", "Mardi"="2", "Mercredi"="3", "Jeudi"="4", "Vendredi"="5", "Samedi"="6", "Dimanche"="7"))
df_ID3$Day <- str_replace_all(df_ID3$Day, c("Lundi"="1", "Mardi"="2", "Mercredi"="3", "Jeudi"="4", "Vendredi"="5", "Samedi"="6", "Dimanche"="7"))
# Heures
df_ID1$Hour <- format(df_ID1$START_TIME, format = "%H")
df_ID2$Hour <- format(df_ID2$START_TIME, format = "%H.%M") #on ajoute les minutes puisqu'il y a 2 obs par heure pour l'ID2
df_ID3$Hour <- format(df_ID3$START_TIME, format = "%H")
# On garde dans une variable juste la date sans l'heure
df_ID1$Date <- as.Date(df_ID1$START_TIME, format = "%m-%d-%Y")
df_ID2$Date <- as.Date(df_ID2$START_TIME, format = "%m-%d-%Y")
df_ID3$Date <- as.Date(df_ID3$START_TIME, format = "%m-%d-%Y")
Les commandes ci-dessus permettent donc de créer les différentes variables. Ce processus d’extraction transforme les variables créées en ‘caractère’, or les modèles GAM ne supportent pas ce format pour les modélisations donc nous les convertissons au format numérique par la fonction as.numeric() seulement pour la variable ‘Hour’ qui prend les heures et les minutes pour l’ID2, sinon nous appliquons la fonction as.integer() pour les autres variables puisque les valeurs sont discrètes.
Nous voilà désormais avec 9 variables par sous-base de chronique. L’intérêt de séparer quelques informations contenues dans la variable ‘START_TIME’ se trouve dans les modélisations GAM : le fait de les avoir passées au format numérique permet d’appliquer des fonctions des variables dans les modèles, et ainsi améliorer les prévisions.
De nombreux packages sous R offrent des statistiques intéressantes et concrètes pour aider à comprendre les données, tels que sumarytools, DataExplorer et explore. Nous regardons les principales informations de notre variable à expliquer grâce à la première librairie, en considérant les temps des jobs par chronique car les écarts sont très importants d’une à l’autre comme nous avons pu le montrer précédemment.
| JOB_D ID1 | JOB_D ID2 | JOB_D ID3 | |
|---|---|---|---|
| Mean | 12.11 | 178.12 | 1973.89 |
| Median | 10.00 | 172.00 | 1825.00 |
| Std.Dev | 7.04 | 22.15 | 437.12 |
| Min | 7.00 | 151.00 | 1100.00 |
| Max | 58.00 | 396.00 | 3622.00 |
Sur les bases propres c’est à dire débarrassées des valeurs atypiques, on constate une fois encore l’intérêt de diviser en 3 bases du fait des écarts importants dans les durées des jobs. Pour pallier à de tels écarts, nous aurions pu standardiser les temps pour être sur la même échelle mais cela aurait compliqué les interprétations et non avons donc choisi de garder 3 bases distinctes.
Avec une semaine prise au hasard parmi les observations de la chronique 1, par défaut avec la fonction ggplot() on voit que l’ajustement est mauvais, d’où l’intérêt de réaliser des modèles GAM où l’on pourra spécifier chaque paramètre de la régression.
Nous allons dans cette partie commencer les modélisations des durées des jobs, en construisant par un modèle simple GAM par chronique puis en essayant d’optimiser les paramètres et ainsi améliorer les prévisions. Pour évaluer la performance des modèles nous regarderons différentes choses :
Aussi, pour comparer les modèles entre eux nous regarderons les critères suivants :
Toutes ces informations auront pour but de constater la qualité de nos prévisions et ainsi sélectionner le meilleur modèle, c’est à dire celui qui maximise le \(R^2\) et minimise l’AIC, le GCV et le MAPE.
Introduit dès les années 80, le modèle additif généralisé (GAM) permet de traiter des relations non linéaires que peuvent prendre les variables explicatives, en s’affranchissant de la linéarité pour mieux prendre en compte les seuils et les valeurs extrêmes. Le principe d’un GAM est ainsi d’expliquer une phénomène non par des variables explicatives mais par des fonctions de ces dernières : appelées fonctions de lissage non paramétriques ou non spécifiques. Effectivement, aucun paramètre n’est nécessaire à la construction de GAM ; la distribution de la variable à expliquer peut alors prendre la forme d’une loi binomiale, d’une loi de Poisson, de Gamma etc. Le modèle linéaire général est ainsi un cas particulier du modèle linéaire généralisé, où Y suit une loi normale et où les fonctions de liaisons avec les variables explicatives sont de simples fonctions identité telles que f(x)=x.
Pour étudier le comportement des fonctions sur Y, il convient de définir et estimer les fonctions non linéaires \(f_{j}\) qui remplacent les coefficients \(\beta_{j}\). Valables pour les variables quantitatives, ces fonctions permettent de maximiser la qualité de la prévision de Y et ainsi augmenter les performances des modèles construits. Dans le cadre de ce travail nous cherchons à obtenir les meilleures prévisions des durées de différents “jobs” en secondes, c’est pourquoi nous aurons recours aux modèles additifs généralisés.
Plusieurs paramètres entrent en compte lors de l’estimation de GAM et leur ajustement peut permettre une meilleure qualité prévisionnelle, nous jouerons donc sur les paramètres suivants :
Par la combinaison de tous ces paramètres on cherche la fonction la plus lisse mais qui soit la plus précise possible pour modéliser Y, on cherche le compromis entre le log de vraisemblance du modèle et sa complexité - compromis définit par le paramètre de lissage \(\lambda\). L’idée d’un GAM est d’utiliser une forme plus flexible de régression entre Y et X ; mais un \(\lambda\) trop petit conduit à une régression trop précise et donc inefficace qu’on appelle over-fiting ou sur-apprentissage (le cas extrême correspondant à une fonction de base par observation). La complexité est aussi observable par l’edf (estimation du nombre de degrés de liberté) qui représente de combien la variable est lissée donc plus l’edf est grand et plus les splines ont une forme complexe.
Pour les modélisations nous procéderons à la même approche pour les 3 ID ; nous modéliserons dans un premier temps un modèle additif généralisé le plus simple possible, en incluant les variables qui vont avec la chronique mais sans préciser de paramètres énumérés ci-dessus. Puis, petit à petit nous effectuerons une recherche du meilleur modèle en testant différents paramètres jusqu’à trouver les prévisions les plus justes.
# On estime un premier modèle gam le plus simple possible c-à-d. sans rien définir
gam1_ID1 <- gam(formula=JOB_DURATION ~ TAP_VERSION + s(Week) + s(Day,k=1) + s(Hour), data = df_ID1, method="REML")
summary(gam1_ID1) #R2=0.348
##
## Family: gaussian
## Link function: identity
##
## Formula:
## JOB_DURATION ~ TAP_VERSION + s(Week) + s(Day, k = 1) + s(Hour)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.46087 0.15400 80.915 <2e-16 ***
## TAP_VERSION1.16.0 -0.23851 0.46818 -0.509 0.6105
## TAP_VERSION1.17.0 0.07533 1.34894 0.056 0.9555
## TAP_VERSION1.17.1 0.15173 1.07576 0.141 0.8878
## TAP_VERSION1.17.2 -0.17043 0.56229 -0.303 0.7618
## TAP_VERSION1.17.3 -1.64044 0.50608 -3.241 0.0012 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 3.350 4.064 6.596 2.21e-05 ***
## s(Day) 1.868 1.982 3.506 0.025 *
## s(Hour) 8.975 9.000 290.297 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.348 Deviance explained = 35.1%
## -REML = 15794 Scale est. = 32.355 n = 4992
Pour ce premier GAM nous incluons les variables créées ainsi que ‘TAP_VERSION’ de la base de la chronique 1 (‘df_ID1’). On voit alors que cette dernière n’est pas très pertinente puisque seule 1 de ses modalités est significative, au seuil de risque de 5%. Par ailleurs, la variable des heures est elle très pertinente et considérer sa relation à Y de manière non linéaire permet effectivement d’améliorer le modèle (p-value de <2e-16 inférieure à 0.05). Le constat est le même pour les fonctions des variables ‘Week’ et ‘Day’ qui sont elles aussi significatives, cependant on voit que le degré de liberté des jours dans le semaine est très proche de 1 signifiant ainsi une relation presque linéaire avec Y.
plot(gam1_ID1,residuals=FALSE, shade=TRUE, pages=1)
L’utilité de considérer des fonctions des variables explicatives peut se confirmer par les plots des résidus où nous observons clairement des seuils et des variations pour la variable ‘Hour’, une relation légèrement incurvée pour les numéros des semaines, et une quasiment linéaire pour le jour de la semaine sur les temps d’exécution des jobs.
# On plote les valeurs prédites et les valeurs observées pour voir la qualité du modèle simple
fitted <- gam1_ID1$fitted.values
real_ID1 <- df_ID1$JOB_DURATION
date_ID1 <- df_ID1[,2]
prev <- cbind(date_ID1, real_ID1, fitted)
prev <- data.frame(prev)
ggplot(data=prev[60:227,])+
geom_line(mapping=aes(y=real_ID1,x= START_TIME,color="valeurs observées")) +
geom_line(mapping=aes(y=fitted,x= START_TIME,color="valeurs prédites")) +
theme_bw() +
scale_color_manual(values = c(
'valeurs prédites' = 'blue',
'valeurs observées' = 'red')) +
labs(color = 'JOB DURATION', x = "Temps", y = "Durée des jobs",
title = "Valeurs prédites vs. valeurs observées du GAM n°1")
Sur le graphique ci-dessus nous pouvons constater les valeurs prédites par le modèle (en bleu) et les valeurs observées de Y (en rouge) ; nous voyons que ce premier GAM a du mal à modéliser une si grande volatilité et capte mal les valeurs importantes que nous voyons tous les jours à 5h41 pour la chronique 1. On cherche donc a améliorer ce modèle en prenant en compte les constats faits lors de son analyse.
Pour le premier modèle additif généralisé nous avons décidé de ne préciser aucun paramètre pour les fonctions des variables explicatives que nous avons incluses. Dans une deuxième estimation nous définissons le nombre de fonctions de base k par le nombre de valeurs que prennent les différentes variables pour avoir une régression la plus précise possible. Ainsi pour la fonction de l’heure nous posons k=24 puisque les valeurs s’étendent de 0 à 23 heures, de la même manière pour la variable ‘Day’ qui prend des valeurs entre 1 et 7 pour les jours de la semaine nous définissons k=7. Pour la variable ‘Week’ qui indique le numéro de la semaine dans l’année ses paramètres étaient plus difficiles à définir puisque le jeu de données de la chronique 1 est à cheval sur 2 années ; les valeurs de la semaine vont donc de 1 à 52 mais elle ne prend que 30 valeurs différentes (40 à 52è semaine de 2019 et 1 à 18è de 2020). Concernant le type de spline appliqué pour les heures, un comportement cyclique est clairement identifiable donc nous avons ajouté l’argument bs=“cc”, et pour capter la saisonnalité des fonctions de ‘Week’ et ‘Day’ nous avons spécifié bs=“ps” qui indique des p-splines après avoir testé différents types et que celui-ci se soit avéré le plus pertinent.
# Deuxième GAM : on précise les paramètres de celles qui sont utiles
gam2_ID1 <- gam(formula=JOB_DURATION ~ TAP_VERSION + s(Week, k=30, bs="ps") + s(Day, k=7, bs="ps") + s(Hour, k=24, bs="cc"), data = df_ID1, family="gaussian")
summary(gam2_ID1)$r.sq #0.562
## [1] 0.5616597
summary(gam2_ID1)$s.table #2 fonctions significatives
## edf Ref.df F p-value
## s(Week) 3.441082 4.174156 9.900653 5.402969e-08
## s(Day) 2.134445 2.555667 4.123971 9.765975e-03
## s(Hour) 21.928767 22.000000 287.964782 0.000000e+00
plot(gam2_ID1,residuals=FALSE, shade=TRUE, pages=1) # relation quasi linéaire pour 'Day' mais pas pour 'Hour' (cf. edf)
Lorsque nous spécifions les paramètres pour chaque fonction de variable dans le modèle, nous voyons que la qualité du modèle croît beaucoup avec un \(R^2\) plus de 20% supérieur à celui du premier GAM et égal à 56.2%. Par la fonction ‘modele$s.table’ nous voyons la significativité des fonctions des \(X_i\) ; elles le sont toutes au seuil de risque de 5%. Nous constatons aussi que l’edf de la variable des jours de la semaine ‘Day’ est légèrement plus élevé que dans notre premier modèle mais reste proche de la parfaite linéarité (edf=1). Regardons l’évolution des valeurs observées et prédites à nouveau par un graphique temporel.
Comme visible sur le graphique ci-dessus, augmenter le nombre de fonctions de base a permis d’améliorer nos prévisions : nous voyons que les valeurs extrêmes sont beaucoup mieux captées par ce nouveau modèle en revanche les petites fluctuations entre les valeurs importantes sont moins bien modélisées et plus lisses. L’idéal serait donc de trouver le modèle qui capte à la fois les valeurs extrèmes et les variations entre celles-ci.
Dans ces troisièmes estimations nous nous intéressons aux effets non-linéaires bivariés où l’on suppose une corrélation entre certaines variables, ces interactions sont importantes dans les séries temporelles à double saisonnalité. Lorsque l’on inclut dans une même fonction différentes variables explicatives, on suppose que leurs effets se recoupent. Dans un premier temps nous utiliserons une fonction simple d’interaction s(\(X_1\),\(X_2\)) pour voir quelles variables peuvent avoir un effet bivarié, puis on cherchera à préciser cet effet. On ajoute donc au modèle n°2 existant les effets bivariés pour les variables ‘Hour’, ‘Day’ et ‘Week’ en essayant les 3 combinaisons possibles.
# GAM 3 : on ajoute une fonction qui prend un effet bivarié
gam3_ID1 <- gam(JOB_DURATION ~ TAP_VERSION + s(Week, k=30, bs="ps") + s(Day, k=7, bs="ps") + s(Hour, k=24, bs="cc") + s(Hour, Day), data = df_ID1, family = gaussian)
summary(gam3_ID1)$r.sq #0.5616
## [1] 0.5615901
summary(gam3_ID1)$s.table # s(Hour, Day) non significatif
## edf Ref.df F p-value
## s(Week) 3.441641 4.174854 9.9012332 5.390910e-08
## s(Day) 2.134411 2.555628 4.1240830 9.768169e-03
## s(Hour) 21.928036 22.000000 264.7467137 0.000000e+00
## s(Hour,Day) 1.000015 1.000030 0.2120357 6.452041e-01
gam3bis_ID1 <- gam(JOB_DURATION ~ TAP_VERSION + s(Week, k=30, bs="ps") + s(Day, k=7, bs="ps") + s(Hour, k=24, bs="cc") + s(Week, Day), data = df_ID1, family = gaussian)
summary(gam3bis_ID1)$r.sq #0.5617
## [1] 0.5616596
summary(gam3bis_ID1)$s.table # s(Week, Day) non significatif
## edf Ref.df F p-value
## s(Week) 3.441080e+00 4.174153 9.900070e+00 5.408902e-08
## s(Day) 2.134446e+00 2.555668 4.123968e+00 9.766006e-03
## s(Hour) 2.192841e+01 22.000000 2.879624e+02 0.000000e+00
## s(Week,Day) 1.049490e-05 27.000000 2.885171e-07 4.073018e-01
gam3ter_ID1 <- gam(JOB_DURATION ~ TAP_VERSION + s(Week, k=30, bs="ps") + s(Day, k=7, bs="ps") + s(Hour, k=24, bs="cc") + s(Week, Hour), data = df_ID1, family = gaussian)
summary(gam3ter_ID1)$r.sq #0.5651
## [1] 0.5651265
summary(gam3ter_ID1)$s.table # s(Week) non significatif mais s(Week, Hour) significatif
## edf Ref.df F p-value
## s(Week) 1.000001 1.000001 0.2090903 6.475018e-01
## s(Day) 2.118505 2.537928 4.1278985 9.848551e-03
## s(Hour) 21.928301 22.000000 222.5610214 0.000000e+00
## s(Week,Hour) 17.121452 21.307252 2.8290703 1.428118e-05
Les combinaisons s(Hour, Day) et s(Week, Day) ne sont pas concluantes car non significatives, mais l’effet bivarié des variables ‘Week’ et ‘Hour’ est lui pertinent même si la fonction de la variable ‘Week’ seule n’est elle pas significative. Nous voyons aussi que la qualité du modèle ne s’est pas améliorée par rapport au modèle 2 puisque nous n’avons pu définir aucun paramètre pour l’effet bivarié. Le meilleur compromis entre ces 3 modèles estimé est donc le dernier pour lequel nous observons ci-dessous son ajustement aux valeurs réelles de Y.
Le constat est le même que pour le modèle 2 ; les faibles volatilités de la durée des JOBS sont mal captées par le modèle. Ainsi, l’ajout de l’effet bivarié n’a pas eu d’impact significatif sur la prédiction comme on peut le constater à travers les \(R^2\) et le graphique ci-dessus qui sont très proches pour les modèles 2 et 3.
Par rapport à une simple fonction s(), les tensors products (te) ont l’avantage de pénaliser les variables dans des directions différentes pour chaque terme, ce qui permet de traiter des variables qui ont des mesures différentes. On peut aussi spécifier les paramètres pour chaque variable (le nombre et le type de splines). De plus, étant donné que la fonction d’interaction inclut automatiquement les fonctions des variables séparément, il n’est pas nécessaire de les rajouter individuellement dans la spécification du modèle.
## [1] 0.8620482
## edf Ref.df F p-value
## s(Day) 3.542003 5.0000 8.626999 3.654141e-10
## te(Week,Hour) 407.030753 443.9237 70.046469 0.000000e+00
Le modèle incluant à la fois les effets bivariés et univarié est la spécification qui explique le mieux le temps des jobs avec une qualité d’ajustement très élevée atteignant 86.2%. De plus, les variables sont toutes significatives ce qui n’était pas le cas des modèles précédents. Le package mgcv propose plusieurs fonctions d’interactions (te, ti et t2) et nous avons fait le choix de tester la fonction ‘t2’ en plus du ‘te’ qui permet de pénaliser avec des fonctions de pénalisation différentes, mais les résultats étaient quasiment les mêmes que pour ce dernier. Sur le plot on voit que les pics sont bien captés par le modèle, comme les fluctuations entre 2 journées et la légère hausse après le pic donc à 5h41 est elle aussi bien modélisée.
Bien qu’il soit évident que le quatrième modèle généralisé soit le meilleur en termes de qualité prévisionnelle (visuellement et statistiquement: \(R^2\)), nous allons tout de même le vérifier au travers des indicateurs de comparaison présentés précédemment.
| GAM.1 | GAM.2 | GAM.3 | GAM.4 | |
|---|---|---|---|---|
| Nombre d’arguments | 4 | 4 | 5 | 3 |
| Nombre d’arguments significatifs | 3 | 4 | 4 | 3 |
| R2 | 0.3448 | 0.5617 | 0.5651 | 0.862 |
| GCV | / | 21.9 | 21.79 | 7.53 |
| AIC | 31544.84 | 29575.98 | 29550.89 | 24169.54 |
| MAPE | 0.2175 | 0.1483 | 0.1508 | 0.1 |
| Validation | X | ✓ | X | ✓ |
La table ci-dessus nous montre que sur les 4 modèles 2 ont tous les termes significatifs, le modèle n°4 est celui qui minimise les erreurs (GCV, AIC et MAPE) en comparaison à tous les autres ; c’est donc de loin le meilleur modèle pour l’estimation des jobs de la chronique 1. Étant donné que les GAM doivent valider certaines hypothèses sur les résidus, nous allons les vérifier dans la section suivante pour le modèle gardé.
# Résidus meilleur modèle
par(mfrow=c(2,2))
gam.check(gam4_ID1)
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 20 iterations.
## The RMS GCV score gradient at convergence was 2.52454e-05 .
## The Hessian was positive definite.
## Model rank = 701 / 701
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## s(Day) 6.00 3.54 1.06 1.00
## te(Week,Hour) 689.00 407.03 0.93 0.04 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
concurvity(gam4_ID1)
## para s(Day) te(Week,Hour)
## worst 0.9780102 0.6595535 5.912765e+28
## observed 0.9780102 0.3857351 2.750745e-03
## estimate 0.9780102 0.3729913 2.039991e-05
# Plot 3D des valeurs de Y prédites par le GAM final
# on remodélise en ne prenant que 2 arguments pour faire un graphique 3D
gam_Vis <- gam(JOB_DURATION ~ te(Hour, Week, k=c(24,30), bs = c("cc", "ps")), data = df_ID1, family = gaussian)
# visualisation 3D
vis.gam(gam_Vis, n.grid = 50, theta = 35, phi = 32, zlab = "",
ticktype = "detailed", color = "heat", main = "Plot 3D du modèle GAM N°4")
Afin de valider les modèles GAM, les résidus doivent être indépendants et identiquement distribués (iid). Dans le cas des séries temporelles, c’est une hypothèse très forte sur les résidus qui n’est bien évidemment pas respectée en général. Les valeurs des séries temporelles actuelles sont fortement corrélées aux valeurs passées, de sorte que les erreurs du modèle seront également corrélées, c’est ce que l’on appelle le phénomène d’autocorrélation. Cela implique que les coefficients de régression et les résidus estimés d’un modèle peuvent être biaisés. Pour résoudre ce problème, nous pouvons le faire en incluant un processus autorégressif (AR) pour capter l’autocorrélation des erreurs de notre modèle. Nous allons utiliser la fonction gamm() du package mgcv pour réaliser cette nouvelle modélisation du GAM.4 mais cette fois-ci en prenant en compte le phénomène d’autocorrélation.
gam4_ar0 <- gamm(JOB_DURATION ~ TAP_VERSION + s(Day, k=7, bs="cp") + te(Week, Hour, k=c(30,24), bs = c("ps", "cc")),
data = df_ID1,
method = "REML")
gam4_ar1 <- gamm(JOB_DURATION ~ TAP_VERSION + s(Day, k=7, bs="cp") + te(Week, Hour, k=c(30,24), bs = c("ps", "cc")),
data = df_ID1,
correlation = corARMA(form = ~ 1|Week, p = 1),
method = "REML")
Nous pouvons comparer les modélisations, sans et avec prise en compte de l’autocorrélation des résidus avec une ANOVA qui sélectionne le meilleur modèle sur la base du log de vraisemblance généralisé. Les résultats du test ci-dessous nous confirment qu’il était judicieux d’estimer en prenant en compte une possible autocorrélation des résidus de la variable ‘Week’. Cela se voit grâce au BIC inferieur pour le gam4_ar1 ainsi que la p-value significative dans la comparaison des 2 modèles.
## Model df AIC BIC logLik Test L.Ratio p-value
## gam4_ar0$lme 1 11 25382.95 25454.60 -12680.47
## gam4_ar1$lme 2 12 25363.18 25441.35 -12669.59 1 vs 2 21.76292 <.0001
Afin d’observer l’autocorrélation partielle des résidus normalisés il est possible de réaliser le plot suivant :
Ces graphiques nous donnent les corrélations partielles des résidus avec ses propres valeurs laguées. Les valeurs optimales du pACF se trouvent normalement en dessous des lignes en pointillé ce qui n’est pas vraiment le cas ici, même pour le modèle gam4_ar1. Nous allons donc tenter d’améliorer cela.
Pour choisir les paramètres optimaux du nombre de retards à inclure dans le processus autorégressif sur les résidus, il est possible de passer par la fonction auto.arima du package forecast. Cette fonction choisira le meilleur nombre de retards sur la base de l’AIC. Le code suivant permet de faire cela.
## ar1 ma1 ma2 ma3
## 0.92607777 -0.99713982 0.11062114 -0.02357678
D’après la fonction auto.arima, il serait plus pertinent d’appliquer un ARMA(1,3) sur les résidus pour résoudre au mieux l’autocorrélation. Nous devrions vérifier cela en réestimant le modèle et comparant les autocorrélations partielles entre gam4_ar1 et gam4_arma13, mais étant donné le temps conséquent d’exécution, nous avons décidé de présenter les codes sans les exécuter.
#gam4_arma13 <- gamm(JOB_DURATION ~ TAP_VERSION + s(Day, k=7, bs="cp") + te(Week, Hour, k=c(30,24), bs = c("ps", "cc")),
# data = df_ID1,
# correlation = corARMA(form = ~ 1|Week, p = 1, q = 3),
# method = "fREML")
#anova(gam4_ar1$lme, gam4_arma13$lme)
#layout(matrix(1:2, ncol = 2))
#pacf(resid(gam4_ar1$lme, type = "normalized"), lag.max = 48, main = "pACF of gam4 with AR(1)")
#pacf(resid(gam4_arma13$lme, type = "normalized"), lag.max = 48, main = "pACF of gam4 with ARMA(1,3)")
Nous démarrons les modélisations pour la chronique n°2 avec un modèle GAM basique où nous laissons le modèle ajuster ses paramètres seul, c’est à dire avec la méthode “REML”. Pour rappel les connections ID égales à 2 sont observées l’été 2019, la variable ‘TAP_VERSION’ prend seulement 2 valeurs (‘1.4.33’ pour 2366 observations et ‘1.4.34’ pour les restantes). Pour cette sous-base de la chronique n°2 la valeur des heures considère aussi les minutes que nous avons codées comme étant les décimales après la virgule pour pouvoir passer la variable en quantitative précédemment. Dans un premier modèle généralisé nous incluons donc la variable ‘TAP_VERSION’ bien qu’elle ne prenne que 2 valeurs nous supposons qu’elle ne sera pas significative mais cherchons à en être sûrs pour ne pas passer à côté d’une éventuelle information. Nous ajoutons aussi des fonctions des variables ‘Week’, ‘Day’ et ‘Hour’.
##
## Family: gaussian
## Link function: identity
##
## Formula:
## JOB_DURATION ~ TAP_VERSION + s(Week, k = 1) + s(Day, k = 1) +
## s(Hour)
##
## Parametric coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 178.204 0.419 425.309 <2e-16 ***
## TAP_VERSION1.4.34 -1.626 2.014 -0.807 0.42
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Approximate significance of smooth terms:
## edf Ref.df F p-value
## s(Week) 1.010 1.020 36.42 1.53e-09 ***
## s(Day) 1.951 1.998 20.89 9.87e-10 ***
## s(Hour) 8.721 8.977 45.47 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## R-sq.(adj) = 0.164 Deviance explained = 16.8%
## -REML = 11072 Scale est. = 410.5 n = 2497
On voit d’emblée que la variable ‘TAP_VERSION’ n’est pas significative comme nous l’avons supposé précédemment, mais toutes les fonctions incluses elles le sont. De plus, on remarque que l’évolution des variables de semaine et de jours est inverse par rapport aux estimations de l’ID1. La variable ‘Week’ qui était incurvée est désormais parfaitement linéaire (mais on observe tout de même que les semaines fin août ont un temps d’éxécution plus long que les semaines début juillet avec une hausse linéaire) et celle des jours ‘Day’ est maintenant légèrement concave. Leurs degrés de liberté estimés confirment cela car ils sont proches de 1, nous garderons donc pour les estimations suivantes que les variables ‘Day’ et ‘Hour’. Concernant cette dernière on voit que les temps des jobs ne suivent pas les mêmes évolutions dans une journée que l’ID1, ce que nous avions pu noter dès la visualisation des données dans une partie antérieure, et l’on voit que le pic est désormais autour de 9-10h du matin et que l’après-midi les temps augmentent à nouveau pour diminuer en fin de journée. Finalement, on voit que la qualité d’ajustement du modèle est vraiment faible (-20 points par rapport aux GAM de la chronique une) ; le \(R^2\) est de 16.4% pour ce premier modèle, nous cherchons à augmenter le plus possible cette qualité et résumerons notre recherche en un tableau récapitulatif.
Visuellement la précision faible de nos prédictions se confirme avec une courbe des valeurs prédites par le modèle qui ne suit pas les évolutions des valeurs réelles et capte très mal le pic de la matinée à 9h09.
|
Effets univariés seulement
|
Effets bivariés
|
||||||
|---|---|---|---|---|---|---|---|
| GAM.1 | GAM.2 | GAM.3.1 | GAM.3.2 | GAM.3.3 | GAM.3.4 | GAM.3.5 | |
| Termes inclus dans le modèle | |||||||
| 1er terme | TAP_VERSION | s(Day,k=7) | te(Day,Hour,k=c(7,48)) | te(Week,Hour,k=c(7,48)) | te(Day,Hour,k=c(7,48) | te(Day,Hour,k=c(7,48)) | te(Day,Hour,k=c(7,48)) |
| 2è terme | s(Week,k=1) | s(Hour,k=48) | te(Week,Hour,k=c(7,48)) | te(Week,Hour,k=c(7,48)) | |||
| 3è terme | s(Day,k=1) | te(Week,Day,k=c(7,7)) | |||||
| 4è terme | s(Hour) | ||||||
| Nombre de termes significatifs | 3/4 | 2/2 | 1/1 | 1/1 | 1/1 | 2/2 | 3/3 |
| R2 | 0.164 | 0.3576 | 0.3779 | 0.3742 | 0.3722 | 0.4204 | 0.4313 |
| GCV | / | 321.82 | 317.5 | 317.82 | 319.8 | 302.04 | 296.87 |
| AIC | 22127.24 | 21504.82 | 21468.33 | 21471.76 | 21486.68 | 21338.68 | 21295.04 |
| MAPE | 0.0622 | 0.0504 | 0.0488 | 0.0502 | 0.049 | 0.0459 | 0.0457 |
| Validation | X | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ |
Le tableau ci-dessus résume les différentes modélisations que nous avons mises en place pour trouver les meilleures prédictions du temps des jobs sur la chronique 2. Le modèle n°2 reprend les variables significatives du premier GAM en spécifiant le nombre de splines appliquées; 7 pour les jours de la semaine et 48 pour les heures car comme précisé avant, chaque heure comporte 2 observations. La qualité grimpe de 16.4% à 35.76% soit presque 20 points de plus, puis dans un troisième GAM nous ajoutons les effets bivariés que nous appliquons pour chaque variable. Il apparaît alors que les tensors products permettent d’améliorer grandement le modèle et le \(R^2\) le plus élevé que nous obtenons est de 43.13% correspondant au modèle “3.5” où nous avons inclu les effets bivariés des jours/heures, semaines/heures et semaines/jours. En plus de maximiser le \(R^2\), ce modèle minimise l’erreur par CV (GCV), le critère Akaike et le MAPE qui est de 4.57% d’erreur absolue moyenne, ce qui est plus faible que pour les modélisations de l’ID 1 qui variaient entre 10 et 15%. Nous analysons à présent le modèle 3.5 sélectionné.
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 13 iterations.
## The RMS GCV score gradient at convergence was 0.0008146698 .
## The Hessian was positive definite.
## Model rank = 654 / 654
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(Day,Hour) 335.0 75.6 1.06 0.90
## te(Week,Hour) 276.0 44.4 0.99 0.42
## te(Week,Day) 42.0 28.2 0.98 0.29
## para te(Day,Hour) te(Week,Hour) te(Week,Day)
## worst 1.631445e-12 1.0000000 1.0000000 1.0000000
## observed 1.631445e-12 0.8466615 0.4007018 0.5302344
## estimate 1.631445e-12 0.1565188 0.1952045 0.2851099
De la même façon que pour l’ID1, nous pouvons observer la pertinence de rajouter à l’estimation des résidus du modèle GAM un processus autorégressif d’ordre 1. Le résultat du test Anova ci-dessous confirme encore une fois la nécéssité de prendre en compte le phénomène d’autocorrélation dans l’estimation de l’ID2 par le modèle GAM 3.5. Cela s’observe à travers un AIC et un BIC plus faibles pour le modèle gam3.5_ar1 et par le fait que cette différence soit significative car la p-value est très proche de 0.
## Model df AIC BIC logLik Test L.Ratio p-value
## gam3.5_ar0$lme 1 11 21555.53 21619.57 -10766.77
## gam3.5_ar1$lme 2 12 21484.09 21553.94 -10730.04 1 vs 2 73.44542 <.0001
Les plots ci-dessous nous montrent l’autocorélation partielle des résidus normalisés:
La prise en compte de l’autocorrélation n’est pas très convaincante dans le graphique de droite. Il faudrait que les traits verticaux ne dépassent pas le seuil en poitillé bleu pour que les résidus soient considérés comme non-autocorrélés. La solution serait alors d’optimiser les paramètres du modèle autorégressif en augmentant l’ordre du comportement AR ou en ajoutant une partie MA au processus. Dans notre analyse, il nous a semblé nécéssaire d’expliquer comment traiter l’autocorrélation, mais n’avons pas optimisé la prise en compte de l’autoccorélation pour chaque chronique car le principe est le même et le temps d’exécution est très long.
Enfin, c’est avec la chronique correspondant à la connection ID 3 que nous terminons l’analyse et la construction de modèles GAM. Comme observé en partie de visualisation des données, en plus des tendances au sein d’une journée nous observons pour cette chronique des tendances sur une semaine puisqu’on ne retrouve pas de pics de temps plus importants dans la journée le samedi. Nous incluons donc toutes les variables et construisons les modèles en les améliorant au fur et à mesure les prédictions, comme pour les estimations des chroniques 1 et 2.
|
Effets univariés seulement
|
Effets bivariés
|
||
|---|---|---|---|
| GAM.1 | GAM.2 | GAM.3.1 | |
| Termes inclus dans le modèle | |||
| 1er terme | TAP_VERSION | TAP_EXIT_STATUS | TAP_EXIT_STATUS |
| 2è terme | TAP_EXIT_STATUS | s(Week,k=10) | te(Week,Day,k=c(10,7)) |
| 3è terme | s(Week) | s(Day,k=7) | te(Hour,Day,k=c(24,7)) |
| 4è terme | s(Day,k=1) | s(Hour,k=24) | te(Hour,Week,k=c(24,10)) |
| 5è terme | s(Hour) | ||
| Nombre de termes significatifs | 4/5 | 4/4 | 4/4 |
| R2 | 0.639 | 0.681 | 0.885 |
| GCV |
|
62305.42 | 26049 |
| AIC | 22149.91 | 21969.55 | 20544.62 |
| MAPE | 0.104 | 0.0969 | 0.0445 |
| Validation | X | ✓ | ✓ |
Comme visible sur la table ci-dessus nous avons estimé 3 modèles pour cette dernière sous-base de notre échantillon d’entraînement. Dans le premier nous incluons toutes les variables potentielles pour expliquer le temps des jobs de CONNCETION_ID=3 : il en ressort que ‘TAP_VERSION’ n’est pas significative donc nous l’excluons des estimations. Nous voyons aussi que contrairement aux ID 1 et 2, la qualité d’ajustement est directement élevée (63.9%) donc les variations doivent être plus régulières pour que le modèle les capte mieux. Nous améliorons cette première estimation dans une deuxième GAM où nous spécifions le nombre de splines ; de la même manière que les estimations précédentes nous fixons k au maximum c’est-à-dire au nombre de valeurs que prend chaque variable (dans la 3è sous-base il y a 10 semaines). Aussi, comme les modélisations précédentes, le fait de spécifier les types de splines n’améliore pas les estimations voire les diminue (très légèrement). En outre, dans un dernier modèle nous ajoutons les effets bivariés avec les 3 combinaisons possibles des variables explicatives et atteigons alors un \(R^2\) de 0.885, soit le meilleur que nous ayons pu trouver jusqu’ici. Regardons plus précisément les caractéristiques de ce modèle n°3 que nous gardons comme modèle final.
##
## Method: GCV Optimizer: magic
## Smoothing parameter selection converged after 18 iterations.
## The RMS GCV score gradient at convergence was 0.001604292 .
## The Hessian was positive definite.
## Model rank = 439 / 439
##
## Basis dimension (k) checking results. Low p-value (k-index<1) may
## indicate that k is too low, especially if edf is close to k'.
##
## k' edf k-index p-value
## te(Week,Day) 69.0 62.1 1.00 0.54
## te(Hour,Day) 161.0 134.6 1.04 0.87
## te(Hour,Week) 207.0 54.1 1.00 0.50
## para te(Week,Day) te(Hour,Day) te(Hour,Week)
## worst 0.3631692 1.841721e+28 0.9993514 0.9987963
## observed 0.3631692 9.773962e-01 0.8631211 0.7027251
## estimate 0.3631692 2.204587e-01 0.1704511 0.1422472
Pour conclure nous pouvons voir à la fois par le diagnostique des réisuds, les différents graphiques et les taux d’erreur que le modèle estimé pour cette troisième sous-base est le plus précis que nous ayons pu construire jusqu’ici.
Encore un fois nous pouvons observer les critères d’information de l’AIC et du BIC plus faibles significativement pour le modèle gam3_ar1, c’est à dire où est prise en compte l’autocorrélation des résidus.
## Model df AIC BIC logLik Test L.Ratio p-value
## gam3_ar0$lme 1 12 20989.41 21053.78 -10482.71
## gam3_ar1$lme 2 13 20726.02 20795.75 -10350.01 1 vs 2 265.3959 <.0001
Les plots ci-dessous nous montrent l’autocorélation partielle des résidus normalisés :
Pour la chronique 3, nous pouvons observer visuellement une amélioration des résidus lors de la prise en compte de l’autocorrélation par un AR d’ordre 1. En effet les autocorrélations partielles des résidus normalisés sont majoritairement comprises entre les lignes en pointillé. Malgré cela, il serait tout de même judicieux d’optimiser l’ordre du processus AR pour une meilleure prise en compte de ce phénomène comme nous l’avons développé dans la chronique 1.
Pour terminer cette partie de recherche et d’optimisation des modèles additifs généralisés, nous regardons les graphiques d’évolution des valeurs observées et prédites par chaque modèle retenu.
Maintenant que nous avons construit les modèles les plus précis en termes de prévisions nous allons soumettre des prédictions pour la compétition kaggle. Pour cela nous devons mettre la base ‘test’ au même format que la base ‘train’ pour pouvoir appliquer les modèles construits à partir de cette dernière. La base test contient 1700 observations et 5 variables comme précisé en introduction : notre objectif est donc, à partir des données disponibles, de prédire le temps d’éxécution des jobs Stitch.
En reprenant les commandes pour la base d’apprentissage, nous séparons la base test en 3 sous-échantillons à partir de la CONNECTION ID, puis nous créons les variables de semaine, jours et heures pour chacune d’entre elles grâce à la variable ‘START_TIME’. Enfin, nous passons les variables créées au format numérique pour pouvoir appliquer librement les modèles additifs généralisés.
Nous commençons par appliquer le quatrième modèle sélectionné pour la sous-base de la chronique 1, grâce à la fonction predict.gam() contenue elle aussi dans le package mgcv sous R. Nous réalisons ainsi les prévisions pour les 3 base de données test en créant chaque fois une colonne ‘Predicted’ dans l’échantillon concernée, puis nous supprimerons toutes les variables autres que celle-ci et l’ID, pour retrouver la forme du fichier sample_solution que nous importerons sur Kaggle pour voir les résultats de notre estimation, après avoir rassemblé de nouveau les 3 sous-bases.
Les lignes de code suivantes ont permis un tel processus :
#------- ID 1 ----------
# On fait les prévisions pour la base test de la chronique 1 avec le modèle sélectionné pour cet ID
#test_ID1$Predicted <- predict.gam(gam4_ID1, test_ID1) #ne fonctionne pas car les versions TAP sont différentes, nous retirons cette variable et commencons une nouvelle fois
# On reestime le modèle 4 en enlevant 'TAP_VERSION'
gam4_ID1 <- gam(JOB_DURATION ~ s(Day, k=7, bs="cp") + te(Week, Hour, k=c(30,24), bs = c("ps", "cc")), data = df_ID1)
# on fait les prévisions
test_ID1$Predicted <- predict.gam(gam4_ID1, test_ID1)
#------- ID 2 ----------
test_ID2$Predicted <- predict.gam(gam3.5_ID2, test_ID2)
#------- ID 3 ----------
test_ID3$Predicted <- predict.gam(gam3_ID3, test_ID3)
# On remet les bases au format du fichier 'sample_solution'
# on retire toutes les autres variables
test_ID1 <- test_ID1[,-c(1:4,6:8)]
test_ID2 <- test_ID2[,-c(1:4,6:8)]
test_ID3 <- test_ID3[,-c(1:4,6:8)]
# on regroupe les 3 bases un 1 df à soumettre sur kaggle
sample_solution <- rbind(test_ID1, test_ID2, test_ID3)
# on exporte ces 2 df au format csv
rio::export(sample_solution, "./sample_solution.csv")
Après soumission, nous avons obtenons un taux d’erreur MAPE de 14.28%. Notre étude est donc terminée et les prévisions du temps d’éxecution des requêtes Stitch sont satisfaisantes.